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SUMMARY 

The aeroelastic stability of a two-dimensional cascade oscillating in supersonic axial flow is analyzed 
in the time domain, the aeroelastic model consists of a single-degree-of-freedom typical section structural 
model for each blade of the cascade and an unsteady two-dimensional cascade aerodynamic model based 
on the Euler equations. The Euler equations are solved using a time-accurate Alternating Direction 
Implicit (ADI) solution scheme. The aeroelastic equations are integrated in time. The effect of inter- 
blade phase angle is included in the aeroelastic analysis by an appropriate choice of initial and boundary 
conditions. 

Flutter predictions are obtained from the time response of a flat plate cascade in single-degree-of- 
freedom pitching motion. The results correlate well with those obtained from a separate frequency domain 
flutter analysis for all values of interblade phase angles considered. Flutter results are then presented for 
cascades having airfoil sections representative of a supersonic throughflow fan. The validity of the time 
integration method for a cascade of airfoils at various interblade phase angles is demonstrated. 
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g/c 
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speed of sound 

airfoil semichord 

damping coefficient in pitching 

steady pressure coefficient = 2(p - V\)/ P\^\ 

airfoil chord 

moment coefficient about elastic axis 

total energy of the fluid per unit volume 

fluxes of mass, momentum, and energy in x,y direction 

gap-to-chord ratio 

imaginary component of c m 
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I mass moment of inertia about elastic axis 

a 

i incidence angle 

J Jacobian of transformation 

K a spring constant for pitching 

k b reduced frequency based on semichord 

M Mach number 

m mass per unit length 

N number of blades in the cascade 

generalized force in pitching 
q vector of flow variables in x,y coordinates 

r ft radius of gyration about elastic axis 

t time 

t/c thickness-to-chord ratio 

u,v x,y velocity components normalized by the speed of sound 

V total velocity 

V nondimensional speed, V = V^/bo^ 

x,y Cartesian coordinate system 

x' distance measured from leading edge to trailing edge 

a pitching (torsion) displacement, positive nose up, see figure 1(b) 

c* 0 steady component of pitching displacement a 

unsteady component of pitching displacement a 
inlet flow angle 

7 ratio of specific heats 

6 stagger angle 

\i mass ratio 

£ chordwise direction of transformed coordinate system 

rj normal direction of transformed coordinate system 

critical damping factor in pitching 
p air density 

o interblade phase angle 

t nondimensional time 

oj frequency of forced oscillation 

u a pitching frequency, = (K^/I^) 1 / 2 


2 


Subscripts and superscripts: 

0 d( )/* 

n A )/dt 2 

O' d( )/dr 

( )' d ! ( )/dr 2 

F evaluated at flutter 

1 evaluated at inlet 


INTRODUCTION 

Recently there has been interest in designing a turbine engine that operates with supersonic flow 
through the rotor and stator blade rows of the fan. The supersonic throughflow (SSTF) fan is expected 
to provide about a 10-percent decrease in specific fuel consumption, and about a 25-percent reduction in 
propulsion weight, which would lead to a 22-percent increase in aircraft range (ref. l). Other advantages 
of the SSTF include fewer fan stages for a given pressure ratio, less inlet cowl and boundary layer bleed 
drag, better inlet pressure recovery, and better matching of bypass ratio variation to flight speed. NASA 
Lewis Research Center has initiated an exploratory program to investigate the feasibility of the SSTF 
fan (refs. 1 to 4). 

One of the aspects to be considered in the safe design of the SSTF fan is possible aeroelastic insta- 
bility such as flutter. However, no experimental data is currently available and only numerical studies 
have been undertaken. Because of the complexity involved in modeling supersonic throughflow in a three- 
dimensional rotor or stator blade row, only two-dimensional analyses have been attempted in the 
literature. 

In these two-dimensional analyses, the rotor or stator is modeled as a rectilinear two-dimensional cas- 
cade of airfoils with supersonic axial flow (refs. 5 to 8); this is also referred to as a two-dimensional cas- 
cade with a supersonic leading edge locus. Linearized equations (see, e.g., ref. 9) applicable for an 
unloaded, cascade of flat plate airfoils in inviscid flow have been used in these analyses. The motion of 
the airfoils in each mode of the cascade is assumed to be simple harmonic with a constant phase angle 
between the adjacent blades. This assumption leads to an eigenvalue problem in the frequency domain 
and the stability of the system is determined from the eigenvalues (ref. 10). 

The above analyses have neglected the effects of airfoil thickness, camber, and steady-state angle of 
attack on supersonic cascade flutter. One way to include these effects in the cascade flutter analysis is by 
using Computational Fluid Dynamics (CFD) methods. Recent advances in the field of CFD with regard 
to the development of more efficient algorithms and increase in computer speeds and memory have made 
possible numerical studies aimed at computing the complex flow field associated with a supersonic axial 
flow cascade. 

There is no published research on the aeroelastic analysis of a supersonic cascade with supersonic 
axial flow using the CFD approach. However, efforts to predict steady aerodynamic characteristics of 
supersonic compressor airfoils using Euler/Navier-Stokes equations have been reported in reference 11. 
Results from a two-dimensional, finite difference Euler cascade code for the analysis of unsteady super- 
sonic axial flow have been reported in reference 12. The unsteady motion resulting from blade vibration 
is included in the analysis by using a deforming grid technique. This code has been validated by 
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comparing results for a flat plate airfoil with corresponding results from linear theory. The code is suit- 
able for use in the aeroelastic analysis of SSTF fans. 

The objective of the present investigation is to develop an aeroelastic analysis for the prediction of 
flutter in a cascade of airfoils by using a time integration method. The time integration method is neces- 
sary when aerodynamic and structural nonlinearities, such as shock motions, large amplitude oscillations, 
nonlinear elastic behavior, etc. are to be modeled. The classical frequency domain analysis is incapable of 
handling such nonlinearities. In addition, the time integration method has the following advantages. The 
response for both steady and unsteady loading can be obtained. The method gives a realistic simulation 
of the motions of both the fluid and the cascade. It may also allow a considerable saving in computational 
effort when analyzing cascade configurations with several degrees of freedom. Lastly, it allows the ani- 
mation of motion of the cascade and the flow field for better understanding of the physical phenomenon. 

The present study is restricted to a two-dimensional cascade model. A typical section representation 
of each of the blades of the cascade is assumed. Only a single-degree-of-freedom pitching motion is con- 
sidered for the analysis. It is to be noted that the stability of a single-degree-of-freedom system, with no 
structural damping, can be assessed purely from the aerodynamic coefficients. However, the method 
being developed here will also be useful for systems with structural damping and with more than one 
degree of freedom. Note that this is the first instance that the time integration method is being applied 
to investigate aeroelastic stability of a cascade of airfoils. 

The time integration method being used for flutter analysis consists of the simultaneous integration of 
the structural and fluid dynamics equations (ref. 13) of each blade in the cascade. The structural 
dynamics equations are written for each blade in the cascade, assuming that each blade is a typical sec- 
tion with pitching motion. The unsteady aerodynamic load terms that appear in the aeroelastic equa- 
tions are calculated by using the unsteady two-dimensional, compressible, Euler cascade solver of 
reference 12. The loads determined by solving the fluid dynamic equations are used in the aeroelastic 
equations to determine the structural displacements. The displacements calculated from the aeroelastic 
equations are then used as input to the fluid dynamics equations to determine the loads. These two steps 
are repeated successively. In general, the displacements show an oscillatory variation with increasing 
time. If the amplitudes of the displacements grow in time, then flutter is said to occur. 

A brief description of the method of analysis, the computational grids, and boundary conditions is 
given in the next section, followed by results and discussion. Time responses (variations of pitching 
amplitude with time) have been calculated for a flat plate cascade and a cascade of SSTF fan airfoils in 
supersonic axial flow. The flutter speeds calculated from the responses are compared with those obtained 
from a conventional frequency domain flutter analysis; the aerodynamic coefficients required for this fre- 
quency domain analysis are obtained by forcing the airfoils in harmonic motion. 

At this time, there is no experimental data available to compare with the present calculations. 


ANALYSIS 

The rotor or stator is modeled as a two-dimensional rectilinear cascade of airfoils. Figure 1(a) shows 
the geometry and nomenclature for a tuned cascade; a tuned cascade is one in which all the blades have 
identical structural properties. 
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Aeroelastic Model 


A typical section with only pitching motion is considered for each blade in the aeroelastic analysis. A 
typical section model of a blade is shown in figure 1(b). The typical section has pitching (a) displace- 
ment. The governing aeroelastic equation, which is applicable to each blade in the cascade, can be 
written as 


l a a + C a a + K a « = Q q ^ 

where I is the mass moment of inertia about the elastic axis, is structural damping parameter, 
is the pitching spring constant, is the generalized force (moment about the pitching axis), and the 
dots over a denote differentiation with respect to time, t. Defining r^ = I ft /mb 2 , uj 2 ^ = K^/I^, 

C = C /2(jjA ^ , where r is the radius of gyration, m is the mass, b is the semichord, is the 
natural pitching frequency, and is the critical damping factor in pitching. Equation (1) can be 
written as 
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The generalized force is related to the moment coefficient, c m , about the pitching axis (elastic axis), 
as follows: 
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( 3 ) 


where V x and p l are the velocity and density evaluated at the inlet. Defining, nondimensional time 
and velocity as r = ta 1 /2b and V = V^/ba where is the speed of sound at the inlet, and sub- 
stituting equation (3) into equation (2), the final aeroelastic equations can be written as 
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where ( )' = d( )/dr, M — Vj/a^ and ji = m/wp^ 2 . 

The aeroelastic equation (eq. (4)) is marched in time using the constant average acceleration method 
(Newmark scheme) (ref. 14). It should be noted that at flutter, reduced frequency based on semichord, 
k F , used in a frequency domain analysis (ref. 10) is related to flutter speed V F and flutter frequency a> F 
as 


kp = 


a; F b 
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( 5 ) 
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Aerodynamic Model 


The moment coefficient (c m ) in equation (4), is calculated by integrating the pressure over the airfoil 
surface for each blade at each time step. The pressure distribution is obtained by solving the unsteady, 
two-dimensional, compressible Euler equations on a body-fitted coordinate system in strong conservation 
form using an alternating direction implicit (ADI) procedure. The present Euler code is based on the 
code developed for isolated airfoils in reference 15. A brief outline is given here. Further details of the 
formulation can be found in reference 12 . 

In the Cartesian coordinate system, the two-dimensional unsteady Euler equations may be written as 

q t + F x + G y = 0 ( 6 a) 
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( 6 b) 


Here u and v are the Cartesian components of fluid velocity; p is the density, p is the pressure, e is 
the total energy of the fluid per unit volume, terms F and G are the flux terms along the x and y 
directions, and 7 is the ratio of specific heats; the subscripts denote differentiation with respect to time, 
and spatial coordinates x and y. 

The body-fitted coordinate system (f, 77 , t) is related to the Cartesian coordinates (x, y, t) (fig. 2 ), 
according to the following one-to-one relationship: 

£ = £(x, y, t) T] = r/(x, y, t) t = ta iy / 2 b (7) 

In the (£, r/, r) coordinate system, the two-dimensional unsteady Euler equations may be written as 

q r + = 0 ( 8a ) 
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where 


q = J" 


The Jacobian of the transformation J is given by 

J = £ x 7y ~ ” 

and the metrics of the transformation are given by the 

£ x = J y»,; £ y = - J v 

Once a grid has been constructed, the metrics of the transformation are evaluated numerically. Standard 
central differences are used to compute the quantities such as x^, y £, etc., and these quantities in turn 
are used in equations (9) and (10) to compute £ x , £ y , etc. 

Since equation (8) is highly nonlinear, a stable and efficient solution procedure is required. In the 
present work, the Beam-Warming algorithm (ref. 16) is used with modifications noted by Sankar and 
Tang (ref. 15). Artificial dissipation is added in both directions of the computational plane to increase 
solution stability. The solution is second-order accurate in space and first order accurate in time. 

Grids . — The flow solver uses a C-grid generated using the GRAPE (GRids about Airfoils using 
Poisson’s Equation) code (ref. 17) with modifications for cascade geometries (ref. 11). Multiple blade 
grids are constructed by stacking the individual C-grids for each blade. The flow solver works with one 
grid at a time. For a cascade geometry, the location of the grid points on the interface boundary must be 
known to pass information to adjacent grids across this common boundary (fig. 2(c)). In order to avoid 
difficult interpolations and to facilitate the application of boundary conditions at the interface, it is desir- 
able that grid points along the outer boundaries of adjacent grids line up exactly. This is achieved by the 
stacking of the individual C-grids generated using the GRAPE code. 

The grids are generated only once around the steady positions of the airfoils. Since the airfoils move 
during the computation (unsteady solution), it is necessary to move the surrounding grid too. A simple 
rigid body rotation of the grid (as is commonly done for isolated airfoils) is unsuitable for a cascade. A 
special procedure is required to keep the outer boundary of the grid fixed while allowing the airfoil to 
move. 

A deforming grid technique (ref. 12) is used in this analysis to locate the position of the grid points as 
a function of time. The inner boundary follows the blade motion, while the outer boundary remains fixed 
in space (figs. 2(a) and 3). All interior grid lines connecting the inner and outer boundaries are allowed 
to deform. The amount of deformation is a function of the distance from the surface of the airfoils 
together with a weighting function. The displacement of a grid point (i j) is defined as: 
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relationships: 


v x = - J yfi v v = Jx / 
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where Ax— and Ay— are the spatial displacements of grid point (ij), and Axj j 
spatial displacements of grid point (ij) over one time step assuming that the entire 
rigid body. The weighting function, Wj j is defined as: 
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where s— is the arc length from the airfoil surface (j = l) to some grid point along a i = constant grid 
line, and’ j = jmax corresponds to the outer boundary grid line. From equations (11) and (12), it can be 
seen that grid points on the inner boundary (s— = 0) give w— = 1, which means the airfoil surface fol- 
lows the rigid body motion of the blade. Conversely, the grid points on the outer boundary give w— = 0 
and the positions of these grid points remain fixed at the initial specified locations. The interior grid 
points shear in space relative to the initial grid as w— varies between 0 and 1. The grid velocities can 
be easily found by dividing the grid deformation by tlie time step value. Using this deforming grid 
approach, the metrics must be recalculated at each time step. 

Boundary conditions . — The present flow solver independently solves the flow equations for each grid 
around a blade and uses interface boundary conditions to model the cascade effects. Ghost points are 
assigned at the first interior grid line (j = jmax - l) and are used by the adjacent grid. Although it is 
tempting to use the most current flow information as it becomes available from the integration scheme, it 
is important to only use flow information from the same time step across the interface boundaries. This 
eliminates possible time inaccuracy in multiple blade solutions due to the particular sequence in which 
individual grid solutions are obtained. The metric data are also forced to be continuous along the inter- 
face boundaries. This procedure essentially makes the blade-to-blade interface boundaries invisible to the 
flow solution. Periodic boundary conditions (explained in the next section) are applied at the lower and 
upper boundaries of the cascade (fig. 2(c)). 

The remaining boundary conditions include those at the inlet and exit planes, the blade surfaces and 
the slits aft of the airfoils (fig. 2). At the inlet the flow variables (density, two velocity components, and 
pressure) are specified to be uniform. At the exit these flow variables are extrapolated from the interior 
by using a simple first order model. These boundary conditions are valid for supersonic flow at the inlet 
and exit. Solid wall boundary conditions are applied along the airfoil surface and the flow variables are 
averaged across the slit aft of the airfoil. 

Periodic boundary conditions and interblade phase angle (<y) . — The flow through a two-dimensional 
rectilinear cascade is a model for the flow through a N bladed rotor or stator of a turbomachine. Hence, 
there exists a fundamental spatial periodicity across the ends (in the direction of the stagger line) of the 
cascade consisting of N blades. Periodic boundary conditions can, therefore, be applied at the upper 
and lower boundaries of the N blade cascade; the periodic boundary conditions constrain the flow 
variables to be equal at corresponding points on the two boundaries. However, under certain circum- 
stances it is possible to apply the periodicity condition across a much smaller domain, as will be explained 
in this section. Clearly, this can result in a substantial saving of computational resources. 
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For steady flow through a stationary cascade, blade-to-blade periodicity of flow variables exists 
naturally. Hence, the periodic boundary conditions are applied at the upper and lower boundaries of a 
single grid; that is, the flow variables are required to have identical values at corresponding points on the 
outer boundary of a single grid (fig. 4(a)). The computations are performed on only one grid rather than 
on N grids since the flow field would be reproduced identically in the remaining grids. 

For unsteady flows in which all blades have identical motions, the blade-to-blade periodicity of steady 
flows is preserved. Hence, the periodic boundary conditions can again be applied at the upper and lower 
boundaries of a single grid . 

For unsteady flows in which all blades have the same periodic motion, except for a constant time-lag 
between adjacent blades, the periodicity condition can be applied across a group containing between 1 
and N grids. Consider a case where all blades move in a simple harmonic motion of the form sin(u;t), 
except for a constant time-lag between adjacent blades. The time-lag can be represented in terms of the 
interblade phase angle ( a ) so that the motion of the N blades is of the form sin(c<;t), sin(o;t + a), 
sin(ci;t + 2a), etc. The interblade phase angle cannot assume arbitrary values, but is restricted to N dis- 
crete values, or = 27rn/N, n = 0, 1, 2,. . ., N - 1, where N is the number of blades of the cascade. This is 
referred to as Lane’s criterion (ref. 18). For the motion described, it is possible to apply the periodicity 
conditions across a group of grids; the number of grids in the group is between 1 and N, depending on 
the value of the interblade phase angle. For example, an interblade phase angle of 180° requires a group 
with two grids (fig. 4(b)), 90° requires four grids (fig. 4(c)), etc. The number of grids in a group (NGG) 
for a given a is given by 


NGG = smallest integer[360/a, 360/(360 - a), N] 
or 

NGG = smallest integer[N/n, N/(N - n), N] 

The use of NGG grids in the computation instead of N grids represents a saving in computational 
effort while providing the correct treatment of boundary conditions for this case of simple harmonic 
motion with constant interblade phase angle. This method of using the required number of grids (NGG) 
for a given interblade phase angle is referred to as the multipassage approach (ref. 19). 

Clearly, a group with the number of grids equal to the number of blades on the rotor can be used to 
compute all the interblade phase angles that arise in flutter calculations (fig. 4(d)). However, the use of 
more than NGG grids for a particular value of a would require additional computational effort with- 
out providing any additional information; the flow field in the group containing NGG grids would be 
reproduced periodically in the remaining (N - NGG) grids. It can also be pointed out that one can use 
more than one group in the computation, so that the total number of grids is between NGG and N. 

For example, for an eight-bladed rotor, and for a = 180°, one can use two, four, six, or eight grids in the 
computation, with periodic boundary conditions applied across the corresponding number of grids. It can 
be verified that the solution in the basic group containing two grids (NGG = 2) will be same in all the 
cases, and this solution will be reproduced in the additional grids being used in the last three cases. 

Periodic boundary conditions and time response calculations . — The preceding discussion was based on 
the assumption of periodic blade motions. In the time integration method, the blade motion is not known 
in advance, and is in general, oscillatory but not periodic. Hence, no fixed frequency or interblade phase 
angle can be identified. The term “interblade phase angle” that will be used in connection with the time 
integration method refers to the interblade phase angle of the forced simple harmonic motion that is used 
to provide the appropriate initial conditions to the time integration method. 
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The time integration calculations would normally be performed over N grids. For a tuned cascade, 
for which all blades are identical, the governing aeroelastic equation for each blade (eq. (1)) is the same. 
Hence, the response of the blades can differ only if the initial conditions are different. For example, 
o — 0° (for the forced motion) will result in all blades having identical responses. In this case it is suffi- 
cient to calculate the response of only one blade. Similarly, for a = 180°, it is sufficient to calculate the 
response of only two blades (using two grids), since the response of all the remaining blades will be the 
same as that of either one of these blades. It can therefore be seen that the concept of identifying a group 
containing NGG grids, across which periodicity can be imposed, remains valid for the time integration 
method. This also means that for a given interblade phase angle, only NGG aeroelastic equations 
(eq. (4)) need to be solved. 

It has been pointed out in reference 12 that for supersonic axial flow the computational domain of the 
flow field can be reduced by inspecting the shock structure through the cascade. The disturbances 
generated inside the Mach cone of the leading edge do not influence the flow outside the cone, which 
means that only a few blades (which fall within the Mach cone of the reference blade) need to be con- 
sidered in the solution for any interblade phase angle. It was also noted that only the aerodynamic loads 
on the reference blade are valid for the specified interblade phase angle. This procedure works quite well 
when the blade motion is simple harmonic with constant interblade phase angle, since the loads on the 
other blades can be calculated once those on the reference blade are known. However, in time-domain 
aeroelastic analysis, the aerodynamic loads on all blades cannot be related to the loads on the reference 
blade, even for a tuned cascade. In this situation, the computations must be made on a group of blades 
depending on the interblade phase angle. By doing this the forces on all the blades and the aeroelastic 
analysis will be accurate. 

Procedure for aeroelastic analysis . — In general, the aeroelastic equation (eq. (4)) is integrated in time 
for each blade independently. As mentioned earlier, the number of blades used in the calculation depends 
on the interblade phase angle; the minimum number of blades required is NGG. The time integration of 
the aeroelastic equations can be started from prescribed initial conditions, a(0) and a'(0) for each blade 
or from initial forced harmonic motion of three to four cycles. In the present study, the latter is followed. 
First, the aerodynamic solution is obtained for three cycles of forced oscillation for given flow conditions 
and interblade phase angle with an assumed reduced frequency. Then, with the motion and the forces (in 
this case the moment) at the end of the third cycle as the initial conditions, the aeroelastic equations are 
integrated (ref. 20). This procedure has three advantages: 

(1) The results can be checked by comparison with those obtained in reference 12. 

(2) The forces can be Fourier decomposed, and a frequency domain flutter analysis can be performed 
(ref. 21). 

(3) It provides the correct initial values for a(0) and a'(0) that correspond to a specified interblade 
phase angle of forced motion. 

To begin the response calculations, a value of nondimensional speed V (eq. (4)) is assumed. Variation 
of the blade motion with time is plotted and stability is assessed from the response. If the response grows 
with time, the value of V is reduced and vice versa. The aeroelastic equations are then integrated 
again. The process is repeated until a value of V is found for which the response is neutrally stable; 
this value of V is the flutter speed. 

Programming aspects . — Multiblade solutions for oscillating cascades are obtained by generating grids 
for each blade. These grids have the same outer boundary shape but different airfoil positions. The 
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solver generates these grids before the unsteady calculation begins by deforming the mean flow grid for 
one blade through a cycle of forced oscillation and storing the grids that occur at multiples of the desired 
interblade phase angle. For example, when o — 90°, four grids (NGG = 4) need to be assembled with a 
blade to blade phase angle of 90° (fig. 4(c)). With the forced oscillation motion assumed to be given by 
a = a Q + sin(ort), where a Q and a l are the mean and amplitude of oscillation, and w is the fre- 
quency of oscillation, the grid from the mean flow solution is used for the first blade (art = 0 ). The grids 
for the remaining three blades are found by deforming the mean flow grid through one cycle of oscillation 
and storing the grid coordinates when art = 90°, 180°, and 270°. The solver numbers the blades and 
stacks the grids to give one group of grids containing four blades. This group is then used as the initial 
grid for the unsteady solution, where it is deformed for the oscillating cascade using the method described 
in the previous section. 


RESULTS AND DISCUSSION 


The response curves (variations of pitching amplitude with time) are presented for a flat plate cascade 
and for airfoil sections of a supersonic throughflow fan rotor and stator undergoing pitching motion. The 
responses are obtained by integrating the aeroelastic equations (eq. (4)) for each blade. The flutter 
speeds are calculated from the response data. As mentioned earlier, first the blades are oscillated in 
forced motion and then equation (4) is advanced in time. Flutter predictions from frequency domain 
solutions are also presented for comparison. 


The results are presented for the following combinations of parameters (ref. 8 ): As previously 
described, a C-grid is generated for each blade by using GRAPE code. 


( 1 ) flat plate airfoil and rotor airfoil 
gap-to-chord ratio, g/c = 0.311 
stagger angle, 6 = 28° 
mass ratio, — 456 

radius of gyration (mid chord), r^ = 0.4311 
computational grid, 199x22 (£ x 77 ) 


( 2 ) stator airfoil 

gap-to-chord ratio, g/c = 0.295 
stagger angle, 0 = 11.7° 
mass ratio, /z — 456 

radius of gyration (mid chord), r ft = 0.4311 
computational grid, 225x22 (f X rj) 


A frequency domain flutter analysis based on linear theory for supersonic axial flow (ref. 8 ), indicated 
that for the structural parameters (/*, r ft ) and cascade parameters ( 0 , g/c) selected here, the instability 
was essentially a pitching motion instability. This suggests that the analysis developed here with only 
pitching motion is sufficient for the present problem. All the calculations are done on a CRAY Y-MP 
computer. The solution takes about 2.11 xl 0‘ 5 CPU sec per time step per grid point per blade and 
requires about 1.8 MW storage. A typical response calculation with two blades requires about 
400 CPU sec. 


Flat Plate Cascade 

The geometry of the flat plate cascade is derived from the SSTF rotor. However, calculations have 
been performed for arbitrary values of interblade phase angle, not limited to the values given by Lane’s 
criterion referred earlier. Results are presented for 6 = 28° and inlet flow angle (/?) of 28°, giving a 
steady angle of attack (a Q ) of 0 °. A thickness ratio of 0.5 percent is used and the leading and trailing 
edges are rounded to aid in the generation of the C-grid (fig. 5). The steady pressure distribution 
obtained for this flow condition is presented in reference 12 . 
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The variation of the imaginary part of the moment coefficient (Im{c m }) with interblade phase angle 
(a) is presented in figure 6 for two values of reduced frequency (k b ); the steady flow Mach number at the 
inlet section is 2.61. These results are obtained by forcing the flat plate cascade in harmonic motion for a 
given a and reduced frequency k b . The moment coefficient is calculated from the pressure distribution 
and Fourier transformed to get the real and imaginary coefficients. The pitching axis is at 30 percent of 
chord (x'/c = 0.3), which was found to be the best location for causing flutter in reference 8. The 
frequency domain flutter solution is obtained as a part of the time integration procedure as mentioned in 
the “procedure for aeroelastic analysis” section. Note that for a pitching airfoil, the Im{c m } indicates 
the stability: positive values indicate unstable motion and vice versa. Also, note that the reduced 
frequency (k b ), used in the analysis is based on semichord, b. 

From figure 6, it can be seen that for k b = 1.00, the Im{c m } is positive for a between 110° and 
260°; for k b = 1.50, Im{c m } is negative for all values of o. It can be seen that the interblade phase 
angle with the highest imaginary moment coefficient is 180°; i.e., the critical interblade phase angle for 
this cascade is 180°. The value of reduced frequency at which Im{c m } becomes zero (neutrally stable 
motion), for a = 180°, is calculated by linear interpolation as k b = 1.4. 

To calculate the neutrally stable condition in the time domain analysis, the response for a = 180° is 
obtained for two values of nondimensional speed, V (=V j/bo; ) (eq. (4)), one for a low V (converging 
oscillations) and the other for a high V (diverging oscillations). The flutter speed corresponding to neu- 
trally stable oscillations (zero damping) is calculated from these two responses by linear interpolation; 
this result is then compared with the frequency domain solution obtained above. Note that at flutter V 
and k b are related by k b = (1/V )((«;/<«/) as given in equation (5); here v a is the frequency of the 
typical section in pitch. 

The variation of the pitching amplitude with time is shown in figure 7 for o — 180°. A group con- 
sisting of two grids is used in the numerical calculations. The response calculations for two values of V 
for blade 1 (corresponding to grid l]l are shown in figure 7(a). The response for V = 0.6 shows con- 
verging oscillations, where as for V =1.0, it shows diverging oscillations. The calculated flutter speed 
for zero damping is V = 0.762 and (w/a;^) = 0.995. The value of k b calculated from the response and 
equation (5) is 1.311, which is 7 percent lower than the one obtained previously from frequency domain 
solution. 

Figure 7(b) shows the response obtained from calculations from blade 2 (corresponding to grid 2). 
Comparing the responses in figures 7(a) and (b), it can be seen that both blades 1 and 2 show the same 
response differing only by the phase angle of the initial forced motion. This indicates that the effect of 
the interblade phase angle is accurately accounted for while integrating the aeroelastic equations. The 
calculations were repeated with a group containing four grids. Identical responses were obtained (not 
shown) for the alternate blades, indicating that a group consisting of two grids is enough for a = 180°. 

In order to further validate the time integration method, responses were calculated at other interblade 
phase angles and correlated with the frequency domain results. From figure 6 it can be seen that the flat 
plate cascade in pitching motion was stable for a = 0°, 90°, and 300° for k b = 1.0. Since (w/u;^) is 
nearly equal to one for this case, it is necessary only to show that the corresponding responses at o = 0°, 
90°, and 300° are stable for V =1.0. For these inter blade phase angles, groups consisting of one, four, 
and six grids, respectively, are used. Figure 8(a) shows the responses obtained for all the four interblade 
phase angles studied here, namely a — 0°, 90°, 180°, and 300°. The variation of the imaginary 
moment coefficient with a (from fig. 6) is repeated in figure 8(b). It can be seen that for a — 0°, 90°, 
and 300°, the response is stable whereas for a = 180°, it is unstable. This indicates that for a = 0°, 
90°, and 300°, a higher V would be required for unstable oscillations. This is consistent with obser- 
vations in the frequency domain analysis for k b = 1.0 as shown in figure 8(b). 
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From the above calculations, it is concluded that the aeroelastic equations are being correctly inte- 
grated and the effect of interblade phase angle is being accurately accounted for in the time integration 
method. 


SSTF Airfoil Cascade 

Having established that the time integration method is properly implemented for cascade analysis, the 
flutter speeds of two airfoil sections are calculated in both frequency and time domains and the results are 
correlated. One of the airfoils corresponds to the rotor and the other to the stator section near the mid- 
span of a supersonic throughflow fan. Additional design parameters of the SSTF fan can be found in 
reference 3. 

Rotor section . — The SSTF rotor airfoil section is about 5 percent thick, and has a chord of 4.065 in. 
This corresponds to a midspan span location on a compressor blade design at NASA Lewis. The solu- 
tions are obtained for /? = 36° and 0 = 28°, giving a steady angle of attack (a Q ) of 8° (fig. 1(a)). The 
steady, static pressure obtained for this flow condition is presented in figure 9(a) and is compared with 
that obtained in reference 11. Both show a good qualitative trend. However, the present solver predicts 
the shock to be located about 10 percent downstream of the position predicted by reference 11. The dis- 
crepancy may be due to viscous and quasi-three-dimensional effects that are included in the analysis of 
reference 11, in addition to different grid sizes. 

The rotor has 44 blades. Lane’s criterion indicates that the rotor can become unstable in any of the 
possible 58 interblade phase angle modes. Therefore, the flutter analysis should be performed for 58 inter- 
blade phase angles, some of which may require up to 58 passages (grids). Fortunately, for this supersonic 
throughflow cascade, the critical interblade phase angle lies between 180° and 230° (ref. 8); this allows 
the calculations to be performed with a small number of passages (NGG = 2 to 8). 

For the results that follow, a frequency domain flutter analysis was done first. The time domain solu- 
tion was then obtained for selected interblade phase angles and the results were correlated with the fre- 
quency domain results. The variation of imaginary moment coefficient with interblade phase angle is 
shown in figure 9(b) for two values of reduced frequency. The moment is taken about a point near the 
midchord. From figure 9(b), it can be seen that for k b = 0.25, the Im{c m } is positive for a between 
110° and 260°; for k b = 0.5, Im{c m } is negative for all values of o. The critical interblade phase angle 
moves from 180° for k b = 0.25° to 230° for k b = 0.5. At flutter, the critical interblade phase angle 
(for which the Im{c m } is the greatest) is 180° and the calculated value of k b at which Im{c m } 
becomes zero is 0.29. 

Again, as was done for the flat plate cascade, the response for o — 180° is obtained for two values of 
nondimensional speed (V ). A group consisting of two grids is used in the numerical calculations. * The 
responses for two values of V from blade 1 (solid lines) are shown in figure 10. Response for V = 2.8 
shows converging oscillations, whereas for V = 3.0, it shows diverging oscillations. The calculated flut- 
ter speed is V = 2.915 and [<jj/u = 0.877; the value of k b calculated from equation (5) is 0.301, 
which is 3.8 percent higher compared to the value obtained from the frequency domain solution. * 

Figure 10 also shows the response calculation from blade 2 (dotted lines) for the same values of V . It 
can be seen that the responses of blades 1 and 2 have the same amplitude but differ in phase angle. 

Compared to the flat plate case studied in the previous section, the loaded SSTF airfoil becomes 
unstable at a higher flutter speed. Also, the response curves show the influence of steady aerodynamic 
loading on the response. In the case of the flat plate, the steady aerodynamic load is zero, and the aero- 
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elastic response oscillated about a mean of zero. However, for the airfoil section, the steady loading is not 
zero and as soon as the forced motion ends the mean motion changes to a value dictated by the steady 
aerodynamic loading. Nevertheless, after the transients due to steady loading have decayed, both the 
blades maintain the same interblade phase angle. 

Stator section : The stator has 52 blades. The airfoil section of the stator is about 5 percent thick, 
and has a chord of 3.49 in. This corresponds to a midspan location on a compressor blade designed at 
NASA Lewis. The solutions are obtained for /? = 23.3° and 9 = 11.7°, giving a steady angle of attack 
(a 0 ) of 11.6° (fig. 1(b)). 

The variation of imaginary moment coefficient with interblade phase angle is shown in figure 11 for 
two values of reduced frequency at an inlet Mach number of M x = 2.72. The moment is taken about a 
point near the midchord point. It can be seen that for k b = 0.3, the Im{c m } is positive for a between 
110° and 260°; for k b = 0.4, Im{c m } is negative for all values of a . The value of k b , at which Im{c m } 
becomes zero is 0.389. The critical interblade phase angle is 180°. 

As was done for the rotor airfoil section cascade, first the response for a — 180° is obtained for two 
values of nondimensional speed (V ). The responses are shown in figure 12. Again, a group consisting of 
two grids is used in the numerical calculations,. The response for two values of V from blade 1 (solid 
lines) are shown in figure 12. Response for V = 2.0 shows converging oscillations, where as for 

V = 2.5, it shows diverging oscillations. The calculated flutter speed is V = 2.30 and {u/v a ) = 0.902. 
The calculated k b is 0.392, which is within 1 percent of the value obtained from frequency domain solu- 
tion. Figure 12 also shows the response calculations from blade 2 (dotted lines) for the same values of 

V . It can be seen that the responses of blades 1 and 2 show the same amplitude but differ in phase 
angle. In figure 12, it is also seen that the steady loading causes the response to oscillate about a new 
mean position after the forced motion ends. 

Additional time domain response calculations (not shown here) for rotor and stator airfoil sections 
also showed the behavior predicted by frequency domain flutter calculations for other interblade phase 
angles. 


CONCLUDING REMARKS 

A time response procedure is developed for prediction of flutter of a two-dimensional cascade in super- 
sonic axial flow. The analysis is based on a typical section structural model and uses a finite difference 
code based on the Euler equations. The initial condition used to start the time integration procedure is a 
harmonic motion of each blade with an interblade phase angle. 

The predicted flutter speeds correlated well with those obtained from a frequency domain analysis. 
The expected behavior is reproduced by the time integration method for the interblade phase angles con- 
sidered. The calculations for the loaded airfoil sections showed higher flutter speeds compared to the flat 
plate cascade. This shows the importance of including the effects of shape, thickness and loading in the 
analysis. 

An alternative procedure for simultaneously calculating the response in all interblade phase angle 
modes is under consideration; this consists of using general (random) initial conditions. This eliminates 
the need to of calculate responses for each interblade phase angle mode separately, as done in the present 
paper. It may provide substantial saving in computational time and this procedure is being investigated 
at the present time. 
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(a) Cascade geometry. 



Figure 1 .—Cascade geometry and typical section model. 
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(b) C-grid in computational plane. 
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Figure 2.— C-grid in physical and computational planes and grid boundary definitions. 




(a) Undeformed grid. (b) Deformed grid. 

Figure 3— Deforming grid technique with exaggerated motion. 


17 





PERIODIC 
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(d) a = 360°/N; NGG = N. 


Figure 4.— Relation between minimum number of grids per group (NGG) and interblade phase angle (a) for an N bladed cascade. 
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(a) Global grid. 



Figure 5.— Grid for flat plate cascade. 



Figure 6— Imaginary part of moment coefficient about xYc = 0.3 versus interblade phase 
angle; flat plate cascade; inlet Mach number M 1 = 2.61 ; stagger angle 0 = 28°; inlet flow 
angle 3 = 28°; gap-to-chord ratio g/c = 0.31 1 ; pitching displacement a = 0.1 sin (cut). 
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(a) Blade 1. (b) Blade 2. 

Figure 7 —Response of flat plate cascade (parameters same as for Fig. 6; interblade phase angle a = 180°). 
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(a) Response of flat plate cascade for interblade angles of 0°, 90°, 1 80°, and 300° (parameters same as for Fig. 6; blade 1 , V* = 1 .0). 



(b) Imaginary part of moment coefficient about x7c = 0.3 versus a, Kb = 1 .0. 

Figure 8 —Response of flat plate cascade for all four interblade phase angles and variation of imaginary moment coefficient with interblade phase angle (a). 
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(a) Steady pressure coefficient distribution for rotor airfoil 
configuration. 


(b) Imaginary part of moment coefficient versus a (rotor section airfoil; parameters 
same as for Fig. 6, except inlet flow angle p = 36°. 


Figure 9 — Steady pressure coefficient distribution for rotor airfoil configuration and variation of imaginary moment coefficient with interblade phase angle (<j). 




Figure 10— Response of rotor section airfoil (parameters same as for Fig. 9; interblade phase angle a = 180°; blades 1 and 2). 



Figure 1 1— Imaginary part of moment coefficient versus interblade phase angle (stator 
section airfoil; inlet Mach number M 1 = 2.72; stagger angle 0 = 11.7°; inlet flow angle p = 
23.3°; gap-to-chord ratio g/c = 0.297; pitching displacement a = 0.1 sin (a>t). 
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BLADE 1 
BLADE 2 




Figure 12— Response of rotor section airfoil (parameters same as for Fig. 11 ; interblade phase angle a = 1 80° ; blades 1 and 2). 
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